      program Coupling    
c************************************************************************************
c   Calculates the transfer integral between 2 MOs in two different fragments in a complex (or in a dimer) 
c   Different intermediate calculations:
c          1- interaction energy between two MOs, one on each fragment: <MOi(fragment1)|H|MOj(fragment2)> 
c          2- energy of each MO in one fragment in the presence of the other fragment     
c          3- Checks for the normalization of each OM
c          4- overlap between the two MOs         
c
c      Before executing this code,you must have submited 3 Gaussian calculations:
c             a- one for fragment 1, producing the file "nameoffile1.com.pun"  containing the MO coeffs of fragment 1 (keyword "PUNCH(OM) NOSYMM" and  "scf=NOSYMM")
c             b- one for fragment 1, producing the file "nameoffile2.com.pun"  containing the MO coeffs of fragment 2 (keyword "PUNCH(OM) NOSYMM" and  "scf=NOSYMM")
c             c- one for the complex, producing the files "FILE.23" and "FILE.24" (23 and 24 are numbers introduced freely by the user, see below) 
c                containing the overlap and fock matrices (between OAs).
c                At the keyword line you must add "POP=NBOREAD NOSYMM" and "scf=NOSYMM"  
c                At the end of the file .com of the complex, after the coords and one blanc line, add the line "$NBO SAO=w23 FAO=W24 $END" 
c
c          It is suggested to rename the files FILE.23 and FILE.24 in "nameoffile1.FILE.23" and "nameoffile1.FILE.24"
c      The code reads the following files:
c           1- The file "data.in" created by "ls nameoffile1.FILE.23 nameoffile1.FILE.24 nameoffile1.com.pun nameoffile12.com.pun >data.in". 
c              Attention: the order of files in data.in is mandatory (yoy can modify the code if you need to change this order, see lines 56-65). 
c              So you need to check it, and you might be obliged to reorder the files in "data.in".
c           2- the file "dimensionOM-n1n2plusieursOM.in" containing the numbers:
c                  "n1" corresponding to the Number of the AOs in the basis set of fragment 1
c                  "n2" corresponding to the Number of the AOs in the basis set of fragment 2
c                  NumOM1i,NumOM1f = Numbers indicating all MOs of the 1st fragment between initial (i) and final (f) values, and
c                  NumOM2i,NumOM2f = Numbers indicating all MOs of the 2d fragment between initial (i) and final (f) values,
c                  between which the couplings will be calculated.
c 
c      The results are written in the files 'Resuls.out....' and 't-values.out.....' in two different formats, for different subsequent use. Feel free to modify..
c**************************************************************************************
      implicit none         
      real dist, F, Fev, S, S12, F12, E1, E2, t
      integer  NOM1, NOM2, m, k, i, j, n1, n2, nmax, lmax, l, p, q, ts                         
      integer  nOM, NumOM1i, NumOM1f, NumOM2i, NumOM2f                           
      parameter (nmax=10000,lmax=30000000)                           
      double precision Fok(nmax,nmax), FFok(lmax), OM1(nmax),OM2(nmax)        
       double precision  Rec(nmax,nmax), Over(lmax)
       character*5 Car1, Dim2
       character*11 B              
       character*13 C              
       character*65 fichier1,fichier2,fichier3,fichier4
c**********************************************************************    
c**********************************************************************
      open(11, File='param.in', status='old')
      read(11,*) Car1, n1, Car1, n2, Car1, NumOM1i, NumOM1f, 
     $  Car1, NumOM2i, NumOM2f
c      
c                                                                                        
      B='Resuls.out'
      C='t-values.out'
      open(10, File='data.in', status='old')
      read(10,*) fichier1
c        45 is the overlap matrix
      open(45, File=fichier1, status='old')
      read(10,*) fichier2
c        35 is the Fock matrix
      open(35, File=fichier2, status='old')
      read(10,*) fichier3
c        15 is the file containing OM of the 1st fragment
      open(15, File=fichier3, status='old')
      read(10,*) fichier4
c        25 is the file containing OM of the 2d  fragment
      open(25, File=fichier4, status='old')
      Open(150, File=B)
      Open(160, File=C)
c********************************************************************
c     Read the line-matrices Fock and overlap written by Gaussian09: 
c     for one square-matrix of dimensions 2nx2n, gaussian writes a half-matrix in only one line,
c     which contains ((2nx2n-2n)/2)+2n)=2*n*n+n terms.
c*******************************************************************
      l=(((n1+n2)*(n1+n2))+n1+n2)/2 
c     Caution:    the matrices Fock and overlap have a 1st line
c                 with the job name. I have considered that there is a single word.
c                 If there are more than one, the additional words must be errased.                  
      Read(35,*) Car1
      Read(35,*) Car1
      Read(35,*) Car1
      Read(35,*) (FFok(k), k=1,l)                           
      Read(45,*) Car1
      Read(45,*) Car1
      Read(45,*) Car1
      Read(45,*) (Over(k), k=1,l)                           
c      write (6,*) (FFok(k), k=1,l)
c********************************************************************
c     Transformation of the line-matrix Fock (read previously) in square half-matrix
c     of dimensions 2nx2n
c********************************************************************
       k=0
       m=n1+n2
       do i=1,m  
           do j=1,i
                k=k+1
                Fok(i,j)=FFok(k)
           enddo
       enddo
c**************************************************************************
c     Transformation of the line-matrix Over (read previously) in square half-matrix    
c     of dimensions 2nx2n                                    
c***************************************************************************
       k=0
       m=n1+n2
       do i=1,m
          do j=1,i
              k=k+1
              Rec(i,j)=Over(k)
          enddo
       enddo
c**************************************************************************       
c       Fill by symmetry the empty places in the matrix  Fock
c**************************************************************************
       do i=1,m 
           do j=i,m
                 Fok(i,j)=Fok(j,i)
           enddo          
       enddo          
c*******************************************************************
c      Fill by symmetry the empty places in the matrix  Rec 
c********************************************************************
       do i=1,m
          do j=i,m
                 Rec(i,j)=Rec(j,i)
          enddo
       enddo
c*******************************************************************       
c     Here starts a big loop finishing just before the "end"
c     The following calculations are repeated for all the combinaions between   
c     the MO(NumOM1i-NumOM1f) of the first fragment and
c     the MO(NumOM2i-NumOM2f) of the 2d fragment 
c*******************************************************************       
c**********Begining of the big loop*********************
      do NOM1=NumOM1i, NumOM1f
        do NOM2=NumOM2i, NumOM2f

c********************************************************************
c  1- Read the coefs of MO1 and MO2 (corresponding to fragment 1 and
c     fragment 2. The code reads all the MOs, the 1st one is replaced by the 2d one and so on. This is not the best algorythm.....
c********************************************************************
       Rewind 15    
       Read(15,*) Car1
          do i=1,NOM1   
                  Read(15,*) Car1, Car1, Car1, Car1
                  Read(15,1400) (OM1(j), j=1,n1)
          enddo       
          write(150, *),
          write(150, *),'*****************************************'          
           write(150, *),
          write(150,1401), NOM1,NOM2
          write(150, *),
 1400 format(5D15.8)     
 1401 format('Interaction between MO ',I3,x,'of the 1st fragment & MO',
     $ 2X,I3,2X,'of the 2d fragment')    
       Rewind 25
       Read(25,*) Car1
          do i=1,NOM2
                  Read(25,*) Car1, Car1, Car1, Car1
                  Read(25,1450) (OM2(j), j=1,n2)
          enddo    
 1450 format(5D15.8)     
c**********************************************************
c  2- Calculation of the energy of one MO of the 1st  monomer
c**********************************************************
      F=0.0
      Fev=0.0
      S=0.0
         do  i=1, n1
              do  j=1, n1
                 F=F+OM1(i)*OM1(j)*Fok(i,j)
                 S=S+OM1(i)*OM1(j)*Rec(i,j)
              enddo
         enddo
                      E1=F
                      Fev=F*27.21
c 
       write(150, 1010), S,F,Fev
 1010 format('S(1,1)=',1x,F12.8,5x,'E(OM1)=',1x,F12.8,1x,
     $       'ua',2x, '(',1x,F10.6,1x,'eV)')
c*********************************************************
c  3- Calculation of the energy of one MO of the 2d monomer
c*********************************************************
       S=0.0
       F=0.0
       Fev=0.0
         do  i=1, n2                                        
               do  j=1, n2                                     
                  p=n1+i
                  q=n1+j
                  F=F+OM2(i)*OM2(j)*Fok(p,q)                
                  S=S+OM2(i)*OM2(j)*Rec(p,q)                
              enddo                                          
          enddo                                         
                      E2=F
                      Fev=F*27.21
c     
       write(150, 1001), S,F,Fev
 1001  format('S(2,2)=',1x,F12.8,5x,'E(OM2)=',1x,F12.8,1x,
     $       'ua',2x, '(',1x,F10.6,1x,'eV)') 
c***********************************************************************     
c  4- Calculation of the interaction energy between 2 MOs : <OM(1)|H|OM(2)>
c     and their overlap : <OM(1)|OM(2)>
c**********************************************************************
        S=0.0 
        F=0.0 
        Fev=0.0 
         do  i=1, n1
            do  j=1, n2
               p=n1+j
                    F=F+OM1(i)*OM2(j)*Fok(i,p)
                    S=S+OM1(i)*OM2(j)*Rec(i,p)
            enddo
         enddo
                      F12=F
                      S12=S
                      Fev=F*27.21
c
       write(150, 1002), S,F,Fev
 1002  format('S(1,2)=',1x,F12.8,5x,'F(1,2)=',1x,F12.8,1x,
     $       'ua',2x, '(',1x,F10.6,1x,'eV)')
c**********************************************************************     
c  5- Calculation of the transfer integral t(1,2) 
c***********************************************************************
      t=(F12-((E1+E2)*S12/2))/(1-S12*S12)
      t=t*27.21
       write(150, 1003), NOM1, NOM2, t
 1003  format(4x,'t(',I3,',',I3,')=',1x,F12.8,1x,'eV')
c********************************************************************** 
c************End of the calculations of the big loop******************* 
c********************************************************************** 
       write(160, 1004), NOM1, NOM2, t
 1004  format(4x,'t(',I3,',',I3,')=',1x,F12.8,1x,'eV')
            enddo
       enddo
          write(150, *),'*********************************************'       
      end                                           
      

